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ABSTRACT 

A pseudospectral numerical method is used for the solution of the Navier-Stokes and mass 
transport equations for a sphere in a sinusoidally oscillating flow with zero mean velocity. The 
flow is assumed laminar and axisymmetric about the sphere’s polar axis. Oscillating flow results 
were obtained for Reynolds numbers (based on the free-stream oscillatory flow amplitude) 
between 1 and 150, and Strouhal numbers between 1 and 1000. Sherwood numbers were computed 
and their dependency on the flow frequency and amplitude discussed. An assessment of the validity 
of the quasi-steady assumption for mass transfer is based on these results. 

INTRODUCTION 

It is known that under certain conditions the rates of heat and mass transfer are enhanced by the 
oscillation of the surrounding fluid (Richardson, 1967). This phenomenon is useful in pulsed 
combustion (Sabnis and Lyman, 1982; Nag and Mukherjee, 1975), drying (Lockwood, 1980), and 
absorption of high-intensity sound in particle-laden flows (Lyman, 1976; Scott, 1975). Of practical 
interest is the specific range of Reynolds and Strouhal numbers under which increases in heat 
and mass transfer can be expected to occur. 

Previously, solutions of the flow and mass transfer problems for a sphere in oscillatory flow 
have been obtained by perturbation techniques for the case in which the fluid displacement 
amplitude,^, is much smaller than the sphere diameter, d (Burdukov and Nakoryakov, 1965; 
Riley, 1966). In this limit, transport of heat and mass are dominated by acoustic streaming. In the 
opposite limit, A/d > > 1, one may use available analytical or experimental results for steady 
flow and mass transfer, assuming that at each instant a quasi-steady state exists. But in many cases 
of technological importance, A/d is of order one, and there are no analytical results for that 
problem. 



The experimental data for mass transfer from spheres were reviewed fourteen years ago 
by Al-Taweel and Landau (1976), who summarized the numerous Sherwood correlation formulae 
and attempted a more rational classification of them based (in part) on the A/d ratio. Al-Taweel 
and Landau suggested that the Sherwood number be correlated with a relation appropriate for 
acoutic streaming when A/d < < 1, and that a quasi-steady Sherwood number correlation be used 
for A/d > 0.75. There is, of course, no reason to expect an abrupt transition at this particular 
value of A/d, and the choice of 0.75 was somewhat arbitrary and based on limited data. Fur- 
thermore, the data reviewed were for relatively high Reynolds numbers, because it was convenient 
to use spheres with diameters on the order of centimeters. For small particles, however, the 
Reynolds numbers are smaller, and the effect of vibrations on heat and mass transfer at low and 
moderate Reynolds numbers has received little attention. 

The present work was undertaken in order to calculate the rate of mass transfer^from a 
solid sphere to an oscillating flow in the intermediate range of A/d (0.1 < A/d < 2) and to 
investigate the transition from acoustic streaming to quasi-steady transport. This work was 
intended to complement an experimental program at Syracuse University on combustion of and 
mass transfer from solid spheres in oscillating flows (Sabnis and Lyman, 1982; Drummond, 1981). 
The experimental results indicated that an enhancement in burning rate of spherical coal particles 
was a function of A/d alone up to A/d = 0.94, whereas the data for mass transfer from sublimating 
naphthalene spheres, corresponding to A/d < 0.48, were best correlated by a dimensionless 
variable appropriate for acoustic streaming 

According to dimensional analysis, the Sherwood number can be expressed as a function 
of A/d (which is the reciprocal of the Strouhal number), the Reynolds number based on the 
amplitude of the oscillating velocity, and the Schmidt number. Since Schmidt numbers for diffusion 
into gases (the medium of interest here) are of order one, the most important parameters are the 
Strouhal and Reynolds numbers. The objective of the present work is to determine the dependence 
of the Sherwood number on these parameters by numerical solution of the Navier-Stokes and 
mass transport equations for oscillating flow around a sphere. 

FORMULATION OF THE PROBLEM 

The analysis of the unsteady transport problem is simplified by the assumptions that the gas flow 
is laminar, incompressible, and axisymmetric about the sphere’s polar axis^, and that the gas and 
sublimand-laden vapor each behave according to the ideal gas law. The flow field is assumed 


1 Computations were also carried out for heat transfer, but to save space the results will not be presented here. 

2 It has been observed by Taneda(1956) that even when the flow upstream of a sphere is steady, when the Reynolds 
number is greater than 130, the flow is no longer axisymmetrical. 
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unaffected by the sublimand motion. The particle temperature is spatially uniform and inde- 
pendent of time, as also are the fluid properties. The effect of natural convection on the transport 
process is not incorporated in the present problem formulation, but this does not imply the flow 
is assumed spherically symmetric. 

The vorticity transport equation is used to describe the flow field, and only the azimuthal 
component of vorticity is non-zero for an axisymmetric flow. It is convenient in the numerical 
analysis to introduce a modified vorticity 

% = -exp(z)t,sin0 (1) 

where E; is the dimensionless vorticity, related to the dimensional vorticity <; *, sphere radius a and 
free-stream velocity amplitude U\ by 

la/U, (2) 


and z is the contracted radial variable 


z = In (r/a) 


(3) 


When written in terms of these dimensionless and modified variables, the vorticity transport 
equation is (see Lin and Lee, 1973): 


M 2 d\ 

4 


eexp(-3z) 
2 sin 0 


dz\dQ 


2i;COt0 




exp(-2 z)D 2 % (4) 


where: 

M 2 = Re ■ St 

Re - 2U , a/v 
St = 2a<x>/U x 


( 5 ) 

( 6 ) 
(7) 


and 


D 2 


d 2 d d 2 d 

— + — --cot© — 

<?z 2 dz dQ 2 dQ 


Here, ip is the streamfunction, defined as 

ij> = y 9 /U x a* 


( 8 ) 


(9) 


and satisfies the equation 
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D 2 ip = ^exp(2z) 


( 10 ) 


The dimensionless form of the mass transport equation is (see Chuchottaworn et. ah, 1984): 

stdv | exp(-3z) /^ ay ax A _ 2exp(-2z) / a 2 r f ay | a 2 y t ar^ 

2 ax sinG v ao az dz as J Re • Sc \d z 2 dz ae z co * aoj 

where Y is the concentration of the sublimand divided by the value at the surface, Sc = v JD is 
the Schmidt number, and D is the diffusion coefficient. 

Initial and boundary conditions on the solution in the physical domain (0 < z < z„) shown 


in Fig. la 

are: 



x = 0: 

V = ^ = 0 


(12a) 


Y = 0 

0 < z < z. 

(126) 


Y = 1 

z = 0 

(12c) 

0 = 0 , Jt : 

d \ p 

ip = ? = — = 0, 

T 5 dQ 

0 < Z < 2, 

(12 d) 


ip = 0 


(12e) 


Z = (d 2 y/dz 2 ) z . 0 


(12/) 


Y = 1 


( 12g) 

z = z„: 

1 2 
\p = -exp(2z„)sin 

Gsin t 

(12h) 


^ = 0 


(120 


Y - 0 


(12/) 


Note in Eq.(12f) the absence of the term a / a z since v = 0 at z = 0. 

NUMERICAL METHOD 

In the present work the vorticity and mass transport equations are integrated using a combination 
of the pseudo-spectral and finite-difference methods. The implementation of the pseudo-spectral 
method (PSM) in this work involves expressing the angular dependence of the streamfunction, 
vorticity, and species concentration in terms of a truncated Fourier cosine series. Finite differences 
are used to represent the radial and time derivatives. 
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The assumed functions for ip and ^are expanded in Fourier series on the sphere as follows 
(Orszag, 1974): 


( 41(2, 0, t )1 

U(z.e.t)/ 


sin 2 0 ^ 

n - 0 


F n (z, T) 
<7„(z,x) 


cosrc© 


(13) 


These series automatically satisfy the boundary conditions along the polar axis. When the series 
in Eq.(13) are truncated at n = /V - 1 and evaluated at the polar angle 0* , one obtains: 


Vic 


U(2.e*,T)/ 


= sin 0, 


N- 1 
rt-o 


F n {z,x) 

G„(z,x) 


cosn0 fc 


(14) 


The set of coefficients F n (z, t) and G n (z,x ) represent ip and ^ , respectively, in spectral 
space. The transformation between physical space and spectral space is achieved by first selecting 
a set of evenly spaced collocation points in the polar angle 0 (Figure 1), 


0 


* 



0 < k < N - l 


(15) 


where N = Jt/A0. The purpose of the shift by Jt/2 N is to avoid numerical problems 
associated with the coefficient 1 / sin 0 in Eqs.(4) and (11). The coefficients in Eq.(14) are 
determined by the usual inversion of a Fourier series (Gottlieb and Orszag, 1977) 


where 



N b 


N- 1 

- V < 

i v/ 

1 cos nd k 

’n k-0 

u*. 

( sin 2 0 t 

0, 

n < 


2, 

n = 

0 

1, 

n > 

0 / 


(16) 


(17) 


This result is used to compute the necessary derivatives and their products in the vorticity transport 
equation, shown in Eqn.(4). The z - derivatives are approximated by second-order central dif- 
ference expressions. When a first-order forward difference is used for the time derivative in 
Eqn.(4), then the vorticity transport equation can be written in the form 

5 t+4, -? t = aAxA(S T+ax ) + ( 1 - a)AxA(^) (18) 


where A (^) represents the remaining terms in the equation. For a = 0 the solution is explicit, 
when a # 0, implicit. Although test calculations were performed using explicit and implicit 
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methods, an explicit method was used for the results presented in the present work. The reason 
for this was that the value of the weighting factor for convergence changed with time; a factor 
found to be beneficial at the start of a cycle was potentially detrimental at later times. 

The spectral representation of Eq.(10) is given by Orszag (1974): 

d 2 F d F n- l 

G„exp(2z) = — -f-— --F n (l+n)(2 + «)-- £ pF p (19) 

O 7L O Z p-n+2 

p + n, even. 


As before, second-order central difference formulas are used for the z - derivatives. Equation 
(19) is used to establish the interior values of F n at the new time step from the (known) boundary 
values of F n and the interior G n at the new time step. An iterative solution was used to solve the 
resultant matrix equations. No weighting factors were used to accelerate convergence. The 
solution was assumed to have converged when the difference in the maximum value of F n was 
less than 0.005 from one time step to the next (an error on the order of 0.05% to 0.1%); this 
typically required between 5 and 50 iterations. 

Boundary conditions on the sphere surface and in the free stream are enforced with finite 
differences as shown, for example, in Gottlieb and Orszag (1977). Special attention is drawn to 
the representation of the vorticity boundary condition at the wall, which in this work is given by 


SL-o- 


Vi.*" Vo.* 
2Az 


( 20 ) 


Other finite difference approximations for the vorticity boundary condition at the wall were 
considered (Roache, 1977; Drummond, 1985), but convergence difficulties were experienced with 
them. 

The normalized concentration variable is given by 

N- 1 

Y(z,Q,x) = ^ l/ n (z, T)cosn0 t (21) 

n-0 

The derivatives of Y required in Eq.(ll) are obtained from Eq.(21). Since the computational grid 
for Y remains unchanged from Fig. 1, the coefficients W n are computed from 

Q N- 1 

W n = T7T I T(z ! 0,,x)cosn0 t (22) 

' V 0 n * - 0 
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Substitution of Eq.(21) into Eq.(ll) yields the following expression for the mass transport 
equation. Again, second-order central differences are used for thez - derivatives, and a first-order 
forward difference used for the time derivative. Once those substitutions have been made, Eq.(l 1) 
can be written in the form: 

y T+ax -y x = P(Ax)z.(y x+At ) + ( 1 - p) A xi(y T ) (24) 

The conditions on Y at the sphere surface and in the free-stream are: 

Y 0 .k m 1 .0<fc< N- 1 (24a) 

Y M , k = 0,0<k<N - 1 (246) 

The boundary conditions along the polar axis are automatically satisfied by the assumed angular 
dependence of Eq.(21). 

We summarize the hydrodynamic and species transport solution procedure selected for use 
in the present work as follows: 

1. Establish the computational mesh with 9 from Eq.(15). 

2. Establish initial values of rp , ^ , Y 

3. Compute F n and G n using Eq.( 16). 

4. Calculate interior values of ^at x + A x from the vorticity transport equation, Eq.(l8) 

5. Compute new values of ip at x + A x from Poissons equation, Eq.(19) 

6. Compute new boundary values of ^ from the new interior ip k 

1. Compute new Y at x + A x from the mass transport equation. 

8. Repeat steps 3-7 until the solution has converged. 

In the explicit calculation procedure of the present work the value of the dimensionless time step 
was 0.02. The convergence criterion for mass transfer calculations is discussed below. 

RESULTS 

A series of steady flow calculations were made to test the pseudo-spectral code before it was 
applied to oscillating flow. A drag coefficient of 1.928 for Re =40 computed by the PSM on a 31 
x 30 grid with Az = 0.07, A0 = 6°, and z«,= 2.1, compared favorably with the 1.805 value from 
the FDM work of Lin and Lee (1973) with Az = 0.1, A0 = 6°, and z„= 3.0. The surface average 
Sherwood number computed by the PSM was 7.34 for Re =40 and Sc = 2.54, which compares 
reasonably well with the value of 7.65 computed from the correlation of Chuchottaworn et. at. 
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(1984). More detailed comparisons are given in Drummond(1985), and a recent discussion of 
steady flow drag correlations provided in the work of Oliver and Chung (1987). The intent of 
these computations was simply to compare results of the PSM code to those of the FDM for 
comparable grids. 

Streamfunction 

Oscillating flow computations were made for a vibrational Reynolds numbers (Reynolds 
numbers based on the free-stream velocity, Eq.6) ranging between 1 and 150, and for Strouhal 
numbers between 1 and 1000. A typical sequence of the computed instantaneous streamlines is 
shown in Fig. 2. The lack of published data on the streamfunction or vorticity made an exact 
comparison of the results impossible, but the streamlines are qualitatively similar to the results 
of Obakata et. al.(1978) for unsteady flow around a cylinder. 

It appears the pseudo-spectral method of analysis of the present work is quite effective 
since the results of Fig.2 were obtained with a 21 x 15 grid (Az = 0.1, A0 = 12°, 2.0). Several 

computer runs were made to investigate the time dependence of the results on grid size. Typical 
results for unsteady flow cases were as follows: the value of the streamfunction at z = 0.2, 0 = 
90° , and x = 1.0 dropped by 4.6% when was extended from 2.0 to 3.5, and it decreased by 
3.1% when A0 was reduced from 12° to 6°.^ 

Fig.3 illustrates the result obtained by time averaging the streamfunction over one period 
and is a clear indication of the presence of secondary streaming. Since A/d in Fig. 3 is equal to 
1.0 and the Reynolds number equal to 20, we do not obtain a streaming pattern similar to the 
results of Riley(1966) whose analysis was for small-amplitude flows (A/d << 1) at low (A/2 < < 
1) or high (A/2 > > 1) frequencies. Although Fig.3 is for A/2 = 20, the inner streaming layer is 
not confined to a region close to the surface. This may be attributed to the fact that A/d is not 
small. The lines from the center of the sphere to the centers of the closed loops in Fig.3 occur 
at angles of 45° to 53° from the axis of oscillation and radial distances r/a between 2.3 and 2.5. 
Lane(1955) calculated that iorA/d < < 1 the center of the inner vortex was at 54°43’ and r/a = 
1.3, although he cautioned that his streamfunction plots, while accurate in the 0-direction, were 
very poor approximations in the r-coordinate. 


3 As expected, integral quantities are less affected by these grid variations; the surface average Sherwood number, for 
instance, varied on the order of 0.25% for the grid extremes described. Given the emphasis on mass transfer for the 
present study, and the fact that the finer grid was 3.5 times more CPU intensive, the 15x21 grid was used for most of 
computations in the present work. 
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Mass Fraction Contours 

Mass fraction contours are shown in Fig.4 at various times over two periods of the free- 
stream velocity, for Re = 100 and St = 1. Similar mass fraction contours were obtained for Re 
between 10 and 100 and St between 0.5 and 10. 


Surface-Average Sherwood Number 

The surface average Sherwood number is computed from: 


Sh 



sin 0d0 


( 25 ) 


Fig.5 illustrates the time dependence of Sh for A/d = 1 and various Re for x between 0 and 12. 
Because the initial conditions specify Y = 0 everywhere except at z = 0, there are initially large 
gradients of Y at the surface. This is reflected in very high values of Sh for 0 < x < 1. The same is 
true for the case oiA/d = 0.1, shown in Fig.6. The clear difference in the time dependence of Sh 
for A/d =\ and A/d = 0.1 for x > 2.0 suggests a different mechanism governing mass transport. 
The oscillatory nature of the Sherwood number in Fig.5 reinforces the idea that at large A/d mass 
transfer is controlled primarily by the oscillatory flow. Conversely, the profiles of Fig.6 suggest a 
mass transfer mechanism independent of the primary oscillatory flow, namely that of acoustic 
streaming. In both cases the effect of Re seems to be primarily to shift the curves to higher or 
lower Sh as opposed to being responsible for any change in the basic mechanism for mass transfer. 

The transition in the nature of the time dependence of Sh is easier to see when Re is fixed 
and A/d varied. For the case Re= 20 and 0.1 < A/d < 1, Fig.7 clearly shows the transition from 
mass transfer controlled by acoustic streaming (A/d =0.1) to mass transfer controlled by the 
primary oscillatory flow (A/d = 1.0). It is not possible to discern from these curves a value of A/d 
at which the transition occurs. The transition point is best identified from the time-averaged 
results, as discussed below. 

It is necessary to obtain time-asymptotic values of Sh in order to compare the results of 
numerical analysis with predictions of published correlations. Also, since the latter invariably 
refer to time-averaged mass transfer rates, the time-dependent Sherwood was time-averaged over 
a period of the free-stream flow oscillation which was sufficiently long after startup so that the 
starting transient had died out. This procedure was followed regardless of whether the mass 
transfer process was quasi-steady or due to acoustic streaming. Since the difference in the time 
averages of the Sherwood number over the ninth and tenth cycles was typically on the order of 
one percent, the average of the tenth cycle is reported here. 
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Fig.8 presents the time-averaged Sherwood number as a function of A/d for various 
Reynolds numbers from 1 to 150. The break in each curve corresponds to the transition from 
acoustic streaming to quasi-steady mass transfer. This occurs at values of A/d between 0.1 and 
0.25, considerably lower than the value of 0.75 chosen by A-Taweel and Landau (1976). Also, the 
A/d at which this transition takes place depends on Reynolds number: at Re = 1 it happens at 
A/d = 0.1, whereas at the higher values Re = 40, 100 and 150 it occurs at A/d = 0.25, independent 
of Re. 

COMPARISON WITH EARLIER WORK AND DISCUSSION 

The time-averaged Sherwood number decreases with increasing A/d at smal \A/d. This behavior 
is opposite to that predicted theoretically by Burdukov and Nakoryakov (1965) from a 
boundary-layer analysis in the limit A/d -»0 such that M 2 > > 1. In that limit there is a Stokes 
shear layer of thickness 6 = Jv/co = d/ M at the surface, as well as a steady acoustic streaming 
boundary layer of the same thickness. Burdukov and Nakoryakov consider mass transfer due to 
acoustic streaming only and use the steady velocity at the outer edge of the boundary layer to 
calculate the convection terms in the convective diffusion equation. Their result for the average 
Sherwood number over the surface of the sphere can be written 

Sh = Cj (Z/j /'fuoV') = C l s[We^~Sc~-~A7d 

where Cj is a constant of order one^. Although the result is written in terms of A/d in the second 
expression, following Al-Taweel and Landau (1976), this is somewhat misleading, because there 
is actually no dependence on the sphere diameter. The important parameter in this limiting case 
is the Reynolds number of the streaming flow 

Re str = U\/ o)v = Re- A/d 

In a thorough review of oscillatory flows Riley (1967) pointed out that one must carefully 
distinguish between various possible limiting cases as A/d -> 0. The analysis of Burdukov and 
Nakoiyakov (1965) is not very clear on this point but evidently corresponds to taking the limit in 
such a way that Re S ( r remains finite. According to the above expression for Re S ( r , this implies that 
Re -» 00 as A/d 0. The physical relevance of this limit is difficult to imagine, because the 
assumptions of axisymmetric and laminar flow would be invalid at high Re. 

The present computations were intended to indicate the transition from acoustic streaming 
to quasi-steady mass transfer at finit eA/d rather than to investigate the limit A/d -»0. One reason 
for this is the impossibility of resolving the inner boundary layer when M > > 1 without using a 


4 Burdukov and Nakoryakov (1965) give Cj - 1.3, but there appears to be an error in their calculation of the average 
over the surface of the sphere. The correct value is Cj = 1.89. 
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radial grid spacing for which the computations would become prohibitively expensive. Never- 
theless, the computations were carried out down to A/d = 0.001 (St = 1000), but for fixed Re. 
This corresponds to the physically realizable situation in which the frequency is increased while 
the velocity amplitude t//, sphere diameter and viscosity are held constant. But since Re S ( r -» 00 
as A/d -» 0 at fixed Re, this limit is quite different from that considered by Burdukov and 
Nakoryakov, so the computed Sherwood number has a different dependence on A/d. 

Burdukov and Nakoryakov ( 1965) also conducted experiments in which camphor balls were 
exposed to high-intensity, high-frequency acoustic waves. In these experiments A/d was on the 
order of 0.05 and Reynolds numbers were on the order of 500. Good agreement between theory 
and experiment was claimed, but when the correct value C\ - 1.89 is used the agreement is poor. 

Padmanabha, Nair, and Ramachandran (1970) measured mass transfer from vertically 
vibrating naphthalene spheres to air for 0.375 < A/d < 0.875 and 200 < Re < 2000. Their 
recommended correlation is 


Sh = 0.63 Re 0 ' 64 (A/d/ 63 

Scatter in the data at low Reynolds number casts doubt on the specific A/d dependence for Re 
below 300, however. Since the Reynolds number in their experiments and those of Burdukov and 
Nakoryakov (1965) are well beyond the range 1 < Re < 150 of the present computations, it is 
not possible to make a direct comparison of our results with experiment. In the experiments of 
Drummond (1981) it was found necessary to provide a steady as well as an oscillatory component 
of velocity of the air flow to prevent build-up of sublimand in the test section. Since the steady 
flow introduces additional complications, it was decided to limit the computations to the case of 
zero mean velocity. 
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CONCLUSIONS 


Computations of the flow about and mass transfer from a sphere in oscillatory unsteady flow were 
carried out with the pseudo-spectal method for values of the Strouhal number between 1 and 
1000 (A/d from 0.001 to 1). The purpose was to obtain accurate solutions of the full Navier-Stokes 
and convective diffusion equations without the approximations inherent in previous analyses, 
which were valid only in the limit A/d < < 1. Since the velocity of the steady secondary flow 
(acoustic streaming) is typically an order of magnitude less than the primary flow velocity, an 
indication of the accuracy of the present method is the fact that the time-average of the instan- 
taneous streamfunction gives a flow pattern resembling those calculated or experimentally 
observed by others. There are still some discrepancies, however, such as the departure of the ip 
= 0 streamline from 0 = 0 and 90°. In the present work the centers of the inner vortex of the 
streaming pattern for St - 1 occur at an angle of approximately 50 degrees to the axis of osci llation, 
corresponding well with the theoretical results of Lane(1955) iorA/d < < 1. 

To summarize, the effects of the displacement amplitude A, velocity amplitude U], and 
frequency oo of the oscillatory flow on mass transfer from a sphere of diameter d are as follows: 

1. For a fixed A/d, mass transfer increases as the Reynolds number U id/ v of the 

oscillatory flow increases. (This case corresponds to a flow driven by an oscillating 
piston with fixed stroke but variable frequency). 

2. At a fixed Reynolds number, mass transfer decreases as A/d increases (or co 
decreases). Above A/d - 0.25 mass transfer is virtually independent of A/d, at least 
for the range of Reynolds number considered here (1 < Re < 150). 
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